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Abstract 

A  simplified  isotropic  numerical  treatment  for  solving  the  anisotropic  electron  transport  phenomenon  in  PEM  fuel  cells  has  been  proposed.  In 
order  to  maintain  appropriate  lateral  current  distribution,  the  in-plane  electronic  conductivity  in  the  catalyst  and  gas  diffusion  layers  is  utilized,  while 
an  extra  contact  resistance  is  added  between  the  gas  diffusion  layer  (GDL)  and  the  current-collecting  land  to  compensate  the  reduced  through-plane 
electronic  resistance.  This  simplified  method  is  also  applicable  for  solving  the  anisotropic  heat  transfer  phenomenon  in  PEM  fuel  cells,  and  it 
improves  numerical  convergence  and  stability  in  three-dimensional  large-scale  simulations. 
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1.  Introduction 

Significant  progress  in  numerical  modeling  and  simulation 
of  PEM  fuel  cells  has  been  made  in  the  past  several  years,  and 
many  three-dimensional  multi-physics  models,  which  closely 
integrate  transport  phenomena  with  electrochemical  kinetics, 
have  been  developed  [1-4].  Comprehensive  reviews  can  also 
be  found  in  the  open  literature  [5-7]. 

In  recent  studies  of  Meng  and  Wang  [8,9],  the  electron  trans¬ 
port  phenomenon  in  PEM  fuel  cells  has  been  addressed.  It 
was  concluded  that  by  solving  the  electron  transport  equation 
in  the  catalyst  layers,  the  gas  diffusion  layers  (GDL),  and  the 
current-collecting  plates  on  both  anode  and  cathode  sides,  it 
enabled  further  capabilities  in  modeling  PEM  fuel  cells,  includ¬ 
ing  direct  incorporation  of  the  contact  resistances  between  the 
membrane  electrode  assembly  (ME A)  and  GDL,  and  between 
GDL  and  the  current-collecting  land.  However,  in  these  stud¬ 
ies,  the  electron  transport  process  in  the  catalyst  layer  and  GDL 
was  assumed  to  be  isotropic,  but  in  reality,  the  GDL  material  is 
highly  anisotropic;  for  example,  the  GDL  electronic  conductiv¬ 
ity  in  the  through-plane  direction  is  300-500  S  m_1,  while  it  is 
3000-5000  S  m-1  in  the  in-plane  direction,  showing  an  order  of 
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magnitude  difference.  As  examined  by  Meng  and  Wang  [8],  dif¬ 
ferent  values  of  the  electronic  conductivity  could  result  in  large 
difference  in  the  lateral  electronic  resistance,  and  as  a  result, 
very  different  current  density  distributions.  Therefore,  in  PEM 
fuel  cell  modeling,  the  anisotropic  property  of  the  GDL  material 
must  be  taken  into  account  properly. 

In  theory,  there  is  little  difficulty  in  considering  the 
anisotropic  transport  phenomenon,  but  in  numerical  calcula¬ 
tions,  an  order  of  magnitude  difference  in  the  electronic  con¬ 
ductivity  in  the  anisotropic  GDL  material  could  cause  significant 
difficulty  in  terms  of  numerical  convergence  and  stability,  espe¬ 
cially  in  numerical  implementation  of  a  PEM  fuel  cell  model 
into  the  commercial  CLD  packages,  i.e.  Star-CD  and  Eluent, 
through  their  user-coding  capabilities. 

In  this  technical  note,  a  simplified  method  for  solving  the 
anisotropic  electron  transport  phenomenon  in  a  PEM  fuel  cell 
is  proposed  and  evaluated.  Lurther  extension  of  this  method  for 
solving  the  anisotropic  heat  transfer  phenomenon  is  discussed 
as  well. 

2.  Theoretical  analysis  and  result  discussion 

We  will  first  analyze  the  effects  of  anisotropic  electronic  con¬ 
ductivities  on  the  current  distribution  and  cell  performance.  As 
presented  in  Meng  and  Wang  [8] ,  the  equivalent  lateral  electronic 
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resistance  can  be  expressed  as 


where  the  parameters,  w  and  t ,  represent  the  half-channel  width 
and  GDL  thickness,  respectively,  and  crs  the  electronic  conduc¬ 
tivity.  Based  on  Eq.  (1),  the  lateral  resistance  in  a  typical  GDL 
with  a  half  width  of  500  [xm  and  a  thickness  of  300  [xm  can  be 
calculated  as 


(500  [xm)2 

2  x  5000  (Sm-1)  x  300  [xm 


0.83  cm2. 


If  we  calculate  the  lateral  electronic  resistance  using  the 
through-plane  electronic  conductivity  of  500  S  m-1 ,  we  will  get 


(500  [xm)2 

2  x  500  (Sm-1)  x  300  pan 


8.33  mQ  cm2. 


In  a  numerical  model  assuming  isotropic  GDL  properties, 
as  detailed  in  Meng  and  Wang  [8,9],  the  effect  of  the  two  dif¬ 
ferent  lateral  electronic  resistances  on  current  distribution  is 
clearly  illustrated  in  Lig.  1.  In  the  present  calculation,  a  single¬ 
channel  PEM  fuel  cell  has  been  simulated,  and  the  cell  geometry 
and  other  relevant  physicochemical  parameters  are  presented 
in  Meng  and  Wang  [8,9].  The  v-coordinate  is  in  the  through- 
membrane  direction,  y  along-the-channel  direction,  z  the  lateral 
direction.  As  can  be  seen  in  Lig.  1,  current  distributions  in  the 
lateral  direction  are  completely  different  in  the  two  cases.  Varia¬ 
tions  of  the  lateral  currents  in  the  middle  of  the  cell  length  from 
the  two  cases  are  further  compared  in  Lig.  2.  In  this  figure,  the 
lateral  variation  of  the  solid  line  is  dictated  by  oxygen  supply,  as 
it  is  in  the  same  trend  as  that  without  solving  the  electron  trans¬ 
port,  as  shown  in  Meng  and  Wang  [8].  In  contrast,  the  lateral 
variation  of  the  dashed  line  is  mainly  determined  by  the  lateral 
electronic  resistance. 

Based  on  the  numerical  results,  in  order  to  maintain  a  cor¬ 
rect  current  distribution  in  the  lateral  direction,  the  in-plane 
electronic  conductivity  of  the  GDL  material  must  be  used  in 
a  numerical  model  of  PEM  fuel  cell  assuming  isotropic  electron 
transport.  However,  in  this  method,  the  electronic  resistance  in 
the  through-plane  direction  becomes  incorrect.  Next,  we  will 
examine  the  electronic  resistance  in  the  through-plane  direc¬ 
tion.  As  presented  in  Meng  and  Wang  [8],  it  can  be  expressed 
as 

RtP  =  -•  (4) 

Therefore,  a  typical  resistance  in  the  through-plane  direction 
is 


D  300  pm 

“  500  Sm"1 


=  6  m  <2  cm 


2 


In  an  isotropic  treatment  using  the  in-plane  electronic  con¬ 
ductivity  to  maintain  proper  lateral  current  distribution,  as  dis¬ 
cussed  earlier,  the  electronic  resistance  in  the  through-plane 
direction  becomes 
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Fig.  1 .  Current  distribution  in  the  middle  of  the  membrane,  (a)  using  the  through- 
plane  electronic  conductivity  of  500  S  m- 1 ,  and  (b)  using  the  in-plane  electronic 
conductivity  of  5000  S  m-1 . 


Rtp,2  = 


5000  Sm-1 


=  0.6  mO  cm  . 


(6) 
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Lateral  Distance  (mm) 

Fig.  2.  Variation  of  current  density  in  the  lateral  direction  in  the  mid-cell  length. 


This  electronic  resistance  is  thus  reduced  by  5.4  m£2  cm2.  The 
reduced  electronic  resistance  is  defined  as  Rr, 


Since  the  electronic  resistance  in  the  through-plane  direction 
significantly  influences  cell  performance,  the  reduced  resistance 
must  be  compensated  in  the  isotropic  treatment  using  the  in¬ 
plane  electronic  conductivity.  In  anisotropic  calculations,  since 
the  lateral  electronic  resistance  is  very  small,  as  calculated  in 
Eq.  (2),  the  overall  effect  of  the  electronic  resistance  is  predom¬ 
inantly  one-dimensional,  in  the  through-plane  direction.  The 
through-plane  electronic  resistance  will  mainly  affect  cell  per¬ 
formance  and  the  along-channel  current  variation,  not  the  lateral 
current  distribution.  Therefore,  a  simplified  method  can  be  estab¬ 
lished  by  adding  the  reduced  resistance  to  the  contact  resistance. 
In  this  treatment,  we  can  not  only  maintain  appropriate  lateral 
current  distribution  but  also  obtain  correct  cell  performance. 

Because  the  contact  resistance  between  ME  A  and  GDL  varies 
under  the  gas  channel  and  the  current-collecting  land,  it  is 
appropriate  to  include  the  reduced  through-plane  resistance  into 
the  contact  resistance  between  GDL  and  the  current-collecting 
land.  Furthermore,  since  current  density  increases  when  passing 
through  the  GDL/land  interface,  assuming  a  uniform  current 
distribution  in  the  lateral  direction,  an  extra  contact  resistance 
accounting  for  the  reduced  through-plane  resistance  should  be 
corrected  as 


(W-w) 

W 


where  the  parameter  W  is  the  half- width  of  the  current-collecting 
plate. 

In  practical  fuel  cell  simulations,  the  extra  contact  resistance 
should  be  slightly  adjusted,  like  a  number  of  other  model  param¬ 
eters,  based  on  a  set  of  experimental  data,  including  the  overall 
cell  performance  and  the  along-channel  current  variations,  as 
discussed  in  [10].  Afterwards,  it  should  be  fixed  for  all  other 
simulations. 


R  =  0  niQcm2 
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Fig.  3.  Current  distribution  in  the  middle  of  the  membrane,  (a)  without  contact 
resistance,  and  (b)  with  contact  resistance. 
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(b)  Lateral  Distance  (mm) 


Fig.  4.  Variation  of  current  density  in  the  lateral  direction  at  the  inlet  and  outlet 
regions,  (a)  without  contact  resistance,  and  (b)  with  contact  resistance. 

In  order  to  examine  the  effects  of  the  extra  contact  resistance 
on  the  current  distribution  and  cell  performance,  two  cases  are 
calculated  with  and  without  the  contact  resistance.  The  current 
distributions  are  illustrated  in  Fig.  3,  and  the  variations  of  the 
lateral  current  density  are  compared  in  Fig.  4.  With  a  contact 
resistance  of  10  m^  cm2  added  between  GDL  and  the  current¬ 
collecting  land,  the  cell  performance  decreases  to  0.83  A  cm-2, 
from  an  average  current  density  of  1.0  A  cm-2  without  the  con¬ 
tact  resistance.  The  variations  of  the  current  density  in  the  lateral 
direction  are  very  similar,  especially  near  the  inlet  region;  near 
the  outlet  region,  the  current  density  under  the  land,  calculated 
with  the  contact  resistance,  becomes  higher  because  of  the  low 
cell  performance,  which  leaves  higher  oxygen  concentration  in 
this  region. 

Since  we  have  encountered  numerical  convergence  problems 
in  our  attempts  to  directly  solve  anisotropic  electron  transport  in 
a  three-dimensional  PEM  fuel  cell  model,  which  fully  couples 


various  transport  phenomena  and  the  electrochemical  kinetics 
[5,8,9],  it  seems  that  the  present  simplified  isotropic  method  is 
more  appropriate  for  solving  this  anisotropic  phenomenon.  This 
method  improves  numerical  stability  and  convergence  in  large- 
scale  simulations. 

Because  heat  is  transferred  out  of  a  PEM  fuel  cell  mainly  in 
the  form  of  thermal  conduction  through  the  current-collecting 
land  [11],  this  simplified  isotropic  treatment  could  be  applied  to 
solve  the  anisotropic  heat  transport  phenomenon  in  a  PEM  fuel 
cell,  as  well. 

3.  Conclusion 

In  this  technical  note,  a  simplified  numerical  method  for  solv¬ 
ing  the  anisotropic  electron  transport  phenomenon  in  a  PEM  fuel 
cell  is  addressed.  In  order  to  obtain  appropriate  lateral  current 
distribution,  an  isotropic  treatment  using  the  in-plane  electronic 
conductivity  has  been  proposed.  An  extra  contact  resistance  is 
added  between  GDL  and  the  current-collecting  land  to  com¬ 
pensate  the  reduced  through-plane  electronic  resistance.  This 
simplified  isotropic  method  can  also  be  applied  to  solve  the 
anisotropic  heat  transfer  phenomenon  in  a  PEM  fuel  cell,  since 
heat  is  transported  out  of  the  cell  almost  exclusively  in  the  form 
of  thermal  conduction  through  the  current-collecting  land.  This 
simplified  isotropic  method  for  solving  anisotropic  transport 
phenomena  in  a  PEM  fuel  cell  improves  numerical  convergence 
and  stability  in  three-dimensional  large-scale  simulations. 
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